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Abstract 

Monte Carlo methods are used to study the phase transition in NH4C1( S ) from 
the orientationally ordered 5 phase to the orientationally disordered 7 phase. 
An effective pair potential is used to model the interaction between ions. Ther- 
modynamic properties are computed in the canonical and isothermal-isobaric 
ensembles. Each ammonium ion is treated as a rigidly rotating body and the 
lattice is fixed in the low-temperature CsCl geometry. A simple extension of 
the Metropolis Monte Carlo method is used to overcome quasiergodicity in 
the rotational sampling. In the constant- NVT calculations the lattice is held 
rigid; in the constant-iVpT calculations the lattice parameter is allowed to 
fluctuate. In both ensembles the order parameter rapidly falls to zero in the 
range (200 - 250)K, suggesting that the model disorders at a temperature in 
fair agreement with the experimental disordering temperature (243K). Peaks 
in the heat capacity and thermal expansivity curves are also found in the 
same temperature range. 
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I. INTRODUCTION AND EXPERIMENTAL SUMMARY 



The general study of phase transitions has been of continuing interest to condensed 
matter chemistry and physics for many years. Particular efforts have been made to charac- 
terize and understand second-order (or continuous) phase transitions and associated critical 
phenomena ||T|. Good examples of such continuous phase transitions are found in certain 
atomic and molecular solids when they undergo order to disorder transitions M. These 
order-disorder transitions are probed experimentally by observing anomalies in thermody- 
namic properties like heat capacities and various susceptibilities. Some examples of systems 
that are known to exhibit disordering phenomena are NH4C1( S ), /3-brass, and KCN( S ). 

In this study we focus on the order-disorder phase transition in NHtClfo), an interesting 
case that has received considerable previous theoretical and experimental attention ||. A 
helpful review of the experimental data up to 1978 is given by Parsonage and Staveley @. 
As discussed in Reference 2 ammonium chloride has three phases. The two phases generally 
referred to as 5 and 7 share the CsCl lattice type shown in Figure 1 , whereas the a phase has 
a NaCl lattice type. First measured and reported by Simon in 1922 [|3| and subsequently by 
other workers the <5 ^ 7 transition has been used as a textbook example of an order- 

disorder phase transition From experiment it is known that in the low-temperature 5 
phase, most of the NH^" molecules are in the same orientation, ( the ammonium ions orient 
with their hydrogen atoms pointing towards the "same" crystallographic axes), and the 5 
phase is said to be "orient at ionally ordered." In the 7 phase, the orientations with respect 
to the chloride ions have no long range order ||. This model of the 5 — > 7 phase transition 
has been supported by various spectroscopic measurements, including neutron scattering 
experiments and NMR studies [§-0- 

The origin of the orientational order in NH4CI can be understood in terms of the in- 
termolecular forces of the crystal. There are attractive interactions between the hydrogen 
atoms in an ammonium ion and the adjacent cage of chloride ions (see Figure 1). Any 
isolated ammonium ion in a chloride cage tends to align with either of two equivalent sets 



2 



of four chloride ions. In Figure 1 the hydrogen atoms point toward the set of chloride ions 
that has been shaded. In addition to the interactions of ammonium ions with the near- 
neighbor chloride cage, the ammonium ions interact with each other electrostatically, with 
the lowest-order, non-zero contribution to the interaction arising from the octapole-octapole 
interactions characteristic of tetrahedral charge distributions. The octapole-octapole inter- 
action has an asymptotic decay of R~ 7 (R is the distance between the nitrogen centers on 
any two ammonium ions). There is an absolute potential energy minimum if two interacting 
ammonium ions are aligned, and there is a higher energy local minimum if two ammonium 
ions have opposite alignments. We can think of the two potential minima for the ammonium 
ion in a chloride cage as two available states. These two available states can be mapped 
onto a two-state Ising model by identifying one orientation with an up spin and the other 
possible orientation with a down spin fL5f . The order-disorder transition has often been 
interpreted |16| in terms of such a two state Ising model with a positive ferromagnetic cou- 



pling constant between near neighbor spins. We shall investigate the energetics of this Ising 
picture within a simple pair site potential model in Section [TJ of this paper. 

Although the interpretation of the transition in terms of the two-state Ising model has 
been insightful, it is not a perfect description of the transition. At atmospheric pressure, 
a measurement of the heat capacity C p as a function of temperature shows a distinctly 
A-like shape with a peak at 243K J7|. While this A character has been used to support 
the interpretation of this transition in terms of the Ising model, it is known that there is a 



discontinuous change in the specific volume at this temperature p7|-|19[, as well as a latent 
heat of transition [||||. This transition at atmospheric pressure is, at least, weakly first 
order, not withstanding the qualitative similarity the system has to an Ising lattice. The 
first-order character of the transition has led to descriptions of the system in terms of a 



compressible Ising model |l6| , p0| -[23| 



It has also been established that there is hysteresis associated with the 5 — > 7 phase tran- 
sition in measurements of the specific volume [|TT| , |TE| , the spin-lattice relaxation time [TT| , |r2fl , 



and neutron scattering measurements of the [111] Bragg intensity 1 12| . Pressure-dependent 



studies of this system |) JIT, 19] show that the nature of the transition is a function of pressure. 



At high pressures ( greater than 1400 atmospheres), the transition changes from first order 
to continuous. This transformation in behavior is characteristic of a tricritical point |L|. The 
location of the tricritical point on the phase diagram is known to shift to lower pressures 
when the ammonium ions are deuterated pT8|JT5[ ] , implying the importance of quantum ef- 



fects. To date no complete explanation of the microscopic origin of the tricritical point in 
ammonium chloride and its deuterated counterparts is available. 

The purpose of the work presented here is to provide information leading to a microscopic 
description of the order-disorder phenomena in ammonium chloride that is more complete 
than the traditional Ising picture. Our motivation for the study comes from recent investi- 



gations |24| into the analogues of phase transition phenomena in small clusters. Many small 
clusters have thermodynamic properties as a function of temperature that have been inter- 
preted as analogues of bulk phase transition behavior. As examples of first-order transition 
behavior, small clusters of rare gas atoms exhibit heat capacity anomalies ||25|| , and rapid 
increases in diffusivities [fZ3J] that are characteristic of bulk melting. Recently, the study of 
phase transition analogues in small clusters has been extended to cases that are similar to 
second-order transitions. Lopez and Freeman EH] observed heat capacity anomalies in model 



Pd-Ni alloy clusters, that are reminiscent of the order-disorder transitions known to occur 



in some bulk alloy materials like /3-brass [p7 |. Studies of analogues to second-order phase 



transitions in other kinds of materials continue to be of interest. Given this background, it 
is our intention ultimately to investigate the possibility of orientational disordering transi- 
tions in ammonium chloride clusters. To do justice to such a study a detailed simulation of 
the analogous bulk transition using the same model potential is important. Providing bulk 
information for future work on the cluster systems is an important motivation of the work 
we report. 

A review of some previous computational studies of transitions in molecular and ionic 
crystals has been given by Klein [J2^J . There have also been studies related to the present 
work by Smith |2{| and by Hiiller and Kane |3(| who focused on the orientational motion of 



the ammonium ions and by O'Shea who studied related octapolar solids. More pertinent 
to the present discussion is the work of Klein, McDonald and Ozaki [R3] who studied NH 4 Br, 



KCIO4 and Li2S04 using molecular dynamics methods. Klein et al. introduced a suitable 
order parameter to monitor the disordering transition in ammonium bromide. We shall make 
use of the same order parameter in the work we report. 

The contents and organization of the remainder of this paper are as follows. In the 
next section the model potential is discussed along with the basic structural features of the 
ammonium chloride lattice implied by the model. In Section III we discuss the computational 
details including how we resolve quasiergodicity problems in the simulations. We present 
the computed thermodynamic properties in Section IV including a discussion of the chosen 
order parameter. Of the properties given in Section IV, we show that the model predicts 
no peak in the isothermal compressibility at the disordering transition for the simulation 
sample sizes used in the current work. In Section V we identify the lack of peak in the 
compressibility with finite size effects by presenting results of thermodynamic properties 
predicted in a compressible Ising calculation as a function of sample size. We summarize 
our findings and discuss their significance in Section VI. 



II. THEORETICAL MODEL 
A. Model Potential 

In the present study we use classical Monte Carlo methods to compute thermodynamic 
properties for NH4C1( S ) as it undergoes an order to disorder transition. The application of 
classical mechanics to pair models of ionic solids has been explored by Klein and coworkers 
PS| , |52|] , who have used molecular dynamics methods to study dynamical processes in several 
ionic crystals. Their work gives us good reason to expect that useful information about 
the order to disorder transition in ammonium chloride can be obtained using a simple pair 
model potential dominated by Coulombic contributions. We independently evaluate the 
success of this approach by using the model potential to predict such bulk properties as 
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the lattice constant and the barrier to rotation of a single ammonium ion between the two 
local potential minima. We shall confirm that the model potential yields values for these 
observables that are in reasonable agreement with experiment, and we can expect that the 
potential is sufficiently accurate for our purposes. The approach we use to validate the 
effective pair potential has been advocated elsewhere 



The specification of the potential surface is somewhat simplified by the dynamical model 
we assume. In all calculations reported here, we treat the NfLj" molecular ions as rigidly 
rotating bodies, with the center of mass of each NH^ fixed at the center of each cubic cell. 
This dynamical model makes it unnecessary to specify an internal vibrational potential for 
the NHj molecule. We need only set NHj - NFP+, NHj - CI", and CP - CP interactions 
to specify a useful effective interaction potential. Each of these interactions is further de- 
composed into a total of six pair atom-atom interaction terms. Letting r designate the 
coordinates of all interacting particles in the system under consideration, we assume that 
the total interaction potential U(r) is given by 

i j>i 

where is the distance between particle i and particle j, and u<i is a pair potential of the 
form 

r y ) = Aij exp(-ayry) + + (2) 

' ij ' ij ' ij 

Most of the parameters used to define the potential are presented in Table 1. In addition as 
used elsewhere we set qu = -35 and = —.40. We also set qci = —1.00 The 



charges used on the nitrogen and hydrogen atoms of the ammonium ion were confirmed by 
independent ah initio calculations |36| . In Table 1 we have assumed transferability of most 



of the atom-atom interaction parameters from other, similar systems studied previously by 
Klein et al. and Pettit and Rossky [|3~7|| . For the parameters not available from previous 



work, we use standard combination rules, i.e. relations of the form 



A? - v A*4«, (3) 



Cij — V CuCjj, (4) 

and 

= 2^^" ^ 

The specific sources for the parameters used in Eq.@ are given as footnotes in Table 1. 

As in any calculation, to generate the results that follow both in this and subsequent 
sections, a finite representation of the lattice was used. To reduce the importance of edge 
effects, we also used standard minimum image periodic boundary conditions ||33|| . In the 
main (or central) sample cell we included either 8 ion pairs (48 atoms) corresponding to 
a 2x2x2 lattice, or 27 ion pairs (162 atoms) corresponding to a 3x3x3 lattice. To evaluate 
efficiently the sums over the periodic images we also used standard Ewald methods with 
vacuum boundary conditions |33||. We found the Ewald calculations to be converged by 



setting the decay parameter of the error functions to 5.55 in units of the inverse length 



of a side of the central sample cell [|3l| and including 125 reciprocal lattice vectors in the 
reciprocal lattice sum. 



B. Properties of the lattice within the model potential 

A primary indication of the validity of the parameters listed in Table 1 is given by com- 
paring the lattice constant and cohesive energy at OK predicted by the model potential with 
experiment. Using the parameters in Table 1, we have minimized the energy of ammonium 
chloride using the cesium chloride phase. We have assumed all ammonium ions in the lattice 
are oriented in the same direction with respect to the crystallographic axes in the manner 
shown in Figure 1. The number of ion pairs in the central simulation cell was taken to be 27, 
and the Ewald sums have been evaluated as discussed above. The resulting lattice constant 
ao (the distance between the nitrogen atom on an ammonium ion and an adjacent chloride 
ion) is calculated to be 3.789 A with a cohesive energy of -759.7 kJ mol -1 . The agreement 
with the experimental lattice constant (3.868 A) |39[ and the experimental cohesive energy 
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at 298K (-697 kj mol -1 ) |4(J is within acceptable limits. The differences between the cal- 
culated and experimental values are a result of finite size and thermal effects as well as the 
details of the model potential. 

Another indication of the accuracy of the model potential is given in Figure 2 where 
the potential energy U of the lattice with 27 ion pairs in the central lattice is plotted as a 
function of the rotation angle of the central ammonium ion [see Figure 1]. At an angle of 
= (where we set the zero of energy in this figure) all the ammonium ions in the lattice 
are oriented with respect to a set of crytallographically equivalent chloride ions. The central 
ammonium ion is then rotated as in Figure 1 until an angle of = it/2 where the central 
ammonium ion is oriented with respect to the alternate set of chloride ions in the lattice. 
At this angle the ammonium ion finds a local potential minimum higher in energy than the 
absolute minimum at = 0. The ~18 kJ mol -1 barrier to rotation evident in Figure 2 is 
in good agreement with experimental estimates based on NMR and other data ||. This 
agreement is another indication that the potential model can be expected to be at least 
sufficiently accurate to account qualitatively for the order to disorder phenomenology of the 
system. 

In addition to calculating the barrier to rotation for the ammonium ions, we investigated 
the effective range of the ammonium-ammonium interactions in the lattice. As an initial 
geometry we considered a lattice with all ammonium ions oriented in the same direction with 
respect to a set of crytallographically equivalent set of chloride ions except for the central 
ammonium ion in the simulation cell. This central ammonium ion was oriented to the other 
set of chloride ions (i.e. at an angle of = 7r/2 in Figure 2). We then compared the energy 
of the lattice with a single ammonium ion out of orientation, with the energy obtained by 
rotating another ammonium ion in the lattice at a distance R from the central ammonium 
ion so that the central ion and the additional ion were oriented in the same direction. It 
is important to recognize that the energy calculations were performed without any nearest 
neighbor assumptions about the range of the interactions. The change in potential energy 
AU obtained from this process is given in Figure 3 as a function of R. At distances beyond 
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the first ammonium ion cage, AU becomes constant, and as is evident from Figure 3, the 
interaction between ammonium ions can be well represented by a nearest neighbor model. 



C. Implications 

The simple pair model expressed in Eq.@ gives rise to a stable CsCl lattice with the 
ammonium ions oriented in the same direction in its lowest energy OK structure. The ori- 
entational ordering is a consequence of the attractive octapole-octapole interaction between 
the ammonium ions. For ammonium chloride the octapole-octapole coupling is sufficiently 
weak that interactions beyond the near neighbors can be neglected. 

These preliminary indications lend support to representing the system by a simple Ising 
model with a positive ferromagnetic coupling constant. However, as indicated in the Intro- 
duction the simple Ising picture is not of sufficient complexity to explain the phenomenology 
of the order to disorder transition. The transition from first order to second order behavior, 
characteristic of a tricritical point, provides justification for investigating the transition with 
more detail than the Ising picture. In the next section we develop the necessary Monte Carlo 
tools to investigate the system within the model potential of Eq.(0). By performing simula- 
tions in the isothermal-isobaric ensemble, coupling between lattice motions and rotational 
motions of the ammonium ions will be included. 



III. SIMULATION METHOD 



The details of the Metropolis Monte Carlo method both in the canonical and in the 



isothermal-isobaric ensemble been discussed in many references |33|J4^J43|1 . In this section 
we explain some of the details specific to the simulations performed in the current work, 
and the approach we used to insure that the simulations were done in an ergodic fashion. 

Since the Ewald sums constituted the dominant fraction of the computational effort in 
this work, we examined the consequences of decreasing the number of reciprocal lattice 
vectors. Although the absolute value of the cohesive energy was sensitive to including fewer 
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reciprocal lattice vectors, we found the lattice constant and the barrier to rotation of the 
ammonium ions (See Section [11 B| ) to change by ~0.3 percent when only 8 reciprocal lattice 
vectors (rather than the converged 125) were included. Our interest is in the order to 
disorder transition, and the transition can be expected to be sensitive to the barriers and 
not to the absolute value of the cohesive energy of the crystal. We tested this assumption 
by comparing fluctuation quantities (e.g. heat capacities, compressibilities, etc.) sensitive 
to the location of the transition in the canonical ensemble with 8 and 125 reciprocal lattice 
vectors included. We observed no significant changes in the computed properties. Most of 
the Monte Carlo results were determined with the smaller set of reciprocal lattice vectors. 

The calculations in the canonical ensemble were performed using central cell sizes of both 
8 ion pairs and 27 ion pairs. For each central cell size, the lattice parameters were adjusted 
to give the minimum energy for the orientation having all the ammonium ions in the same 
direction. About each ammonium ion in the lattice we constructed a set of orthogonal 
Cartesian axes. Each Monte Carlo point in the canonical simulations consisted of a rotation 
of a randomly chosen ammonium ion about one Cartesian axis randomly chosen MM. The 



maximum allowed rotation angle was adjusted so that about fifty percent of the moves 
were accepted. As is typical in Monte Carlo simulations, this maximum allowed rotation 
angle was a function of temperature and was found by performing short run experiments at 
each temperature. In addition to these normal Metropolis moves, we found it necessary to 
include moves with larger maximum displacements ten percent of the time. The purpose of 
the moves with magnified displacements (magwalking) was to insure that the ammonium 
ions were given the opportunity to overcome the potential barrier separating the two minima 
in the potential surface (See Figure 2). Consequently, the maximum displacement in the 
magwalking moves was taken to be n/2. 

Evidence that the magwalking scheme discussed in the previous paragraph provides an 
ergodic distribution is given in Figure 4 where the average potential energy of the 8 ion 
pair lattice is plotted as a function of temperature. The data used in generating both 
curves of Figure 4 was initiated with configurations having the ammonium ions in random 
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orientations. Quasiergodicity difficulties [p5j , f45|l can be anticipated at low temperatures 
where all the ammonium ions are expected to have the same orientation with respect to 
the chloride ions. The data in the upper curve of Figure 4 was generated using Metropolis 
Monte Carlo methods with a fixed step size for all moves. The unstable behavior at low 
temperatures is evident. The instability is a result of rapid quenching of the initial random 
orientations of the ammonium ions into local, high-energy minima. This disordered trapping 
causes the quasiergodicity in the sampling of the rotations. In the lower curve of Figure 4, 
the magwalking scheme described in the previous paragraph was used. The sensible behavior 
at low temperatures is clear, and the simple magwalking scheme can be expected to solve 
the quasiergodicity problems in this system. 

In the isothermal-isobaric simulations the Monte Carlo sampling is with respect to the 
distribution 

p(r) = A' 1 exp(-!3U(r) - f3pV) (6) 

where U is the system potential energy, (3 = 1/kgT where T is the temperature and ks 
is the Boltzmann constant, p is the pressure, V is the volume and A -1 is a normalization. 
To perform the isothermal-isobaric simulations in addition to the rotational moves used in 
the canonical study, the lattice parameter (and consequently the system volume) was varied. 
From numerical experiments we found that including volume fluctuations about forty percent 
of the time gave good convergence of the computed properties. The volume fluctuations were 
included with Metropolis Monte Carlo moves again with a maximum displacement chosen so 
that about fifty percent of the attempted fluctuations were accepted. Unlike the rotational 
moves, we found no evidence of quasiergodicity in the volume fluctuations. When a mixture 
of maximum step sizes were used in the volume moves, we found no statistically significant 
changes in the computed thermodynamic properties. 

The volume fluctuations included in this work correspond to a single vibrational breath- 
ing mode for the lattice. Of course, the real lattice dynamics in NH 4 C1 is more complex than 
this simple picture. However, at least we have been able to include important contributions 
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« = 77 hyF (9) 



to the coupling between the rotational and vibrational modes of the lattice. 

In the calculations that follow, we have calculated several fluctuation quantities in addi- 
tion to the energy E and enthalpy H = E + pV of the crystal. These fluctuation quantities 
include the constant volume and constant pressure heat capacities 

Cv/k B = 6 — + p[<U 2 >-<U> 2 } (7) 

and 

C p /k B = f3 2 [< H 2 > - < H > 2 ] (8) 
the isobaric coefficient of thermal expansion 

If , 

V \dT J 

= (k B T 2 V)' 1 [<VH > - <VxH>] (10) 
and the isothermal compressibility 

V \dp ) 

= {k B TV)- l [< V 2 > - < V > 2 } (12) 

In Eq.(^) the notation <> represents averages in the canonical ensemble, and in Eqs.(jBD,(|10D 
and ( |12"D the notation <> represents averages in the isothermal-isobaric ensemble. In the 
results we report in the next section, Cy,C p ,a and k were calculated directly from these 
fluctuation expressions. 

In the calculations that follow, error bars were estimated at the double standard deviation 
level by "binning" the data and estimating the variance of the bin averages about the total 
walker average. In the limit of large numbers, this method is known to be appropriate for the 
correlated data generated by a Metropolis walk. However, for finite Monte Carlo samples the 
suitability of the binning parameters should be checked for a given model system to ensure 
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the reliability of the error estimates. We performed these checks by comparing the computed 
error estimates to those obtained through covariance calculations of the error |J!|[|7||. We 
observed the two methods to give similar results, indicating that suitable binning parameters 
were chosen for this study 

IV. RESULTS 

In this section we provide results of the Monte Carlo studies of the properties of the 
NH4CI lattice as a function of temperature. As expected we shall find features in the ther- 
modynamic properties as a function of temperature that can be associated with a transition 
from rotational order to disorder. To monitor the degree of orientational order in the lattice, 
we use an order parameter introduced by Klein et al [j3^] . To define the order parameter 
we place the origin of a set of Cartesian axes on the nitrogen atom of each ammonium ion, 
and we orient the Cartesian axes so that the z-axis is orthogonal to a square face of chloride 
ions in the chloride cage. The x and y axes are then oriented along the edges of the lattice 
as shown by the coordinates displayed in Figure 1. For each ammonium ion we define 

^ = ^E^M (13) 

where x l , is the x-component of the coordinate of a unit vector pointing from the origin 
toward the z'th hydrogen atom on ammonium ion j and the summation runs over the 4 
hydrogen atoms on the ammonium ion. Mj is defined so that Mj = 1 when ammonium ion 
j is oriented exactly to one set of chloride ions, and Mj = — 1 when it is oriented exactly to 
the alternate set of chloride ions (see Figure 1). Of course Mj is a continuous function of the 
coordinates, so that it takes on values of ± 1 only when the hydrogen atoms point exactly 
to a set of chloride ions. This definition of Mj enables a mapping of the orientations of the 
ammonium ions onto a spin variable. When Mj is positive we can think of ammonium ion 
j as having a positive or "up" spin and when Mj is negative, we can think of ammonium 
ion j as having a "down" or negative spin. To clarify the mapping, in Figure 5 we display a 
particular configuration of the ammonium chloride lattice taken from a canonical simulation 
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with 8 ion pairs in the central cell at 300K. We have placed arrows on each nitrogen atom in 
the lattice to carry information about the algebraic sign of Mj. Positive Mj is represented 
by an up arrow and can be interpreted as an up spin. Similarly, negative Mj is represented 
by a down arrow and can be interpreted as a down spin. To simplify the discussion often 
we shall describe the orientations of the ammonium ions in terms of these spin variables. In 
analogy with the Ising model, we can also define the magnetization per site M of the lattice 
by 

1 N 

M = -J2M j (14) 

j'=i 

where the summation on j is over the N ammonium ions in the lattice. Associated with the 
magnetization is a susceptibility per site x defined by 

X = N(< M 2 > — < M > 2 ) (15) 

where the notation <> in Eq.(|l"5|) represents an ensemble average. We use M as the order 
parameter for the orientational order to disorder transition in the system. 

A. Canonical simulations 

In the canonical simulations for each temperature we included 50000 passes without the 
accumulation of data followed by 200000 passes with data accumulation. Each pass consisted 
of cycling through the ammonium ions in the central simulation cell and attempting to rotate 
each ion once. 

We begin by examining the changes in computed thermodynamic properties that accom- 
pany alterations in the size of the central sample cell. In Figure |6] we present the constant 
volume heat capacity Cy per ion pair in units of the Boltzmann constant as a function of 
temperature. The upper curve gives results when the central cell consists of 8 ion pairs and 
the lower curve gives results for a 27 ion pair central cell. By increasing the size of the 
central cell the width of the heat capacity maximum narrows as expected. 
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By examining configuration files it is possible to verify that the maxima in the heat 
capacity seen in Figure |6| are a result of the rotational disordering of the ammonium chloride 
lattice. We can confirm this interpretation by monitoring the magnetization as a function of 
temperature. In Figure [?| the curve connected by diamonds are computed values of < M 2 > 
as a function of temperature for the 27 ion pair lattice. We have chosen to present < M 2 > 
rather than < M > because in any finite lattice < M >= at any finite temperature. 
In a completely oriented configuration as the temperature approaches zero, M 2 is always 
1 whereas M can be either 1 or -1. Then at low temperatures < M 2 > must approach 
unity, and at high temperatures < M 2 > must approach zero. This behavior is evident in 
Figure 0. At the transition temperature of ~210K < M 2 > changes rapidly. This transition 
temperature matches the peak of the maximum in Cy as a function of temperature. Also 
presented in Figure [7| is the result of mapping the magnetization onto a spin model. The 
data for this spin model are plotted in Figure [7] as points represented by triangles. In the 
spin model we set 

1 N 

M s = -^Sn(Mj) (16) 

where sign(Mj) is the algebraic sign of Mj and M s represents the magnetization of the 
system described by the spin variables. 

The difference between the two curves given in Figure |7| clarifies the extent to which the 
decrease in the magnetization can be attributed to librational modes. Since M is a contin- 
uous function of the coordinates of each ammonium ion, the magnetization will decrease as 
a function of temperature even if all the ions remain ordered. In contrast, M s will change 
only if an ammonium ion in a particular configuration changes its "spin;" i.e. overcomes 
the barrier between the two minima in the rotational potential surface. By comparing the 
decay of the order parameter to the decay of the Ising-like spin parameter in Figure [7], we 
observe that the initial decay of the order parameter is a result of the onset of librational 
motions. At the transition temperature, the order parameter rapidly decays to zero because 
of the loss of long range order of the NH4 orientations. 
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Although in a finite system < M s >=0, in an actual Monte Carlo simulation < M s > may 
differ from zero, because both inverted configurations may be reachable only in Metropolis 
walks that are sufficiently long. The length of the walk required to actually calculate a zero 
value for < M s > will grow both with decreasing temperature and increasing system size. 
The effect of the finite walks can be made apparent by examining the susceptibility x as 
a function of temperature. Such a graph is given in Figure [| for the lattice consisting of 
27 ion pairs in the central cell. The maximum of \ occurs at the same temperature as the 
temperature at which < M 2 > changes rapidly. The fluctuations in x a t this temperature 
are also large. Of course in the limit of an infinite system < M >^ 0, and the susceptibility 
we calculate should approach the infinite system result with increasing system size. 

B. Isothermal-isobaric simulations 

Canonical simulations do not have sufficient flexibility to incorporate coupling of the 
rotational modes of the ammonium ions with the vibrational modes of the lattice. To 
obtain preliminary understanding of the effects of such couplings, we have performed Monte 
Carlo simulations of the thermodynamic properties of the system in the isothermal-isobaric 
ensemble. By fixing the pressure of the system, the isothermal-isobaric simulations more 
closely match the experimental situation than the canonical studies. In exchange for this 
closer connection with experimental data, calculations in the isothermal-isobaric ensemble 
are computationally more demanding than calculations in the canonical ensemble. The 
range of system sizes that can be studied while achieving acceptable levels of convergence is 
restricted. 

The simulations in the isothermal-isobaric ensemble were performed on a 27 ion pair 
representation of the ammonium chloride lattice in the central cell. At each temperature 
300000 Monte Carlo passes were performed without data accumulation followed by about one 
million passes with data accumulation. Each pass consisted of sequential attempted rotations 
of each ammonium ion in the same manner as in the canonical simulations. Additionally, 
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volume fluctuations were attempted in forty percent of the passes. Ten percent of the 
rotational moves used a maximum rotation angle of 7r/2 using the magwalking scheme found 
to be successful in the canonical simulations. The remaining ninety per cent of rotational 
moves used a maximum rotational displacement determined so that about fifty percent of 
the attempted moves were accepted. 

In Figure || we present < M 2 > as a function of temperature calculated both using the 
definition of the order parameter given in Eq.(|l^) (the points represented by diamonds) 
and the projection of the order parameter onto spin variables (the points represented by 
triangles). The projected spin order parameter is the same as that used in the canonical 
simulations and discussed previously. The region of rapid change in the order parameter is 
at a temperature of about 250K. The onset of the transition is more clearly seen in Figure 



10] where we plot the susceptibility as a function of temperature. The maximum occurs at 
a temperature of about 250K, a temperature where the fluctuations in x are large. 

We measure the average volume of the crystal in terms of the lattice parameter clq. The 
average volume is of interest because we would expect to observe a rapid volume change near 
any first order transition. In Figure |TT] the lattice parameter in A is plotted as a function 
of temperature. Although no rapid change in the lattice parameter is seen, there is a clear 
change of slope at the transition temperature. 

In Figure |12] we give the constant pressure heat capacity C p per ion pair, the isothermal 
compressibility k and the isobaric coefficient of thermal expansion a as a function of tem- 
perature. The onset of the maxima in C p and a match the onset of the order to disorder 

' ' 11 a maximum 



transition as identified by the order parameter. From basic considerations j|8 
is expected in k at the transition temperature as well. As is evident from Figure [12|, k is 
found to be a monatomic increasing function of the temperature with no apparent peak. To 
test the sensitivity of the behavior of k to the number of reciprocal lattice vectors included in 
the Ewald sums, we increased the set from 8 to 125 reciprocal lattice vectors. No significant 
change in the compressibility as a function of temperature was observed. We believe that 
the lack of peak is a consequence of the finite size of the central simulation cell, and we give 
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evidence for this in the next section. 



V. A COMPRESSIBLE ISING MODEL 

In the previous section we showed that in the isothermal-isobaric ensemble at the dis- 
ordering transition, peaks were observed in the constant pressure heat capacity and in the 
isobaric coefficient of thermal expansion. The location of the peak maxima were in good 
agreement with the temperature at which there were rapid changes in the order parameter. 
The identification of the peak maxima with the disordering phenomena seems well justified. 

In contrast to the properties discussed in the previous paragraph, no maximum was 
observed in the isothermal compressibility. Since the compressibility should have a specific 



heat-like divergence at the transition temperature in the infinite system [PS], the lack of 
observation is a concern. In this section we show that similar behavior is found in finite 
representations of a compressible two-dimensional Ising model on a square lattice, and the 
lack of peak in the isothermal compressibility we observed for NH 4 C1 may be attributable 
to finite size effects. 

The compressible model we use is based on the two-dimensional Hamiltonian 

+ (17) 



where Si is the spin on site i and can take on values of + or - 1, the notation < ij > on the 
sum represents summation over nearest neighbors only, a, e and a are parameters and R is 
the lattice constant. The first term in the Hamiltonian is the standard spin-spin interaction 
in the Ising model with a lattice-parameter dependent coupling constant. We choose the 
independence of the coupling constant to decay with distance like the octapole-octapole 
interaction in ammonium chloride (-R~ 7 ). The second term in the Hamiltonian is in the 
form of a usual Lennard- Jones potential and provides a balance for the volume fluctuations 
so that the lattice will relax to a physical nearest neighbor distance. In the Ising simulations 
that follow we have attempted to choose parameters for Eq. ( |17|) that mimic the ammonium 
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chloride results of the previous section. For the bulk Lennard- Jones a parameter, we take a 
value that can make the equilibrium lattice distance near that of the real crystal. We then 
take a = Re/2 1 / 6 where R e = 7.31 Bohr. We take e to match the ammonium chloride lattice 
energy of 8380K. Since the barrier to changing the orientation of a single ammonium ion in a 
completely ordered ammonium chloride lattice is about 1000K, we take a = 2000-Rg/iV S pj n , 
where -^Vgpin is the number of spins in the primary cell. 

We determined the properties of the two-dimensional compressible Ising lattice defined 
by Eq.(|l7D using Monte Carlo simulations in the isothermal-isobaric ensemble. The sim- 
ulations consisted of two-million moves without data accumulation followed by 10 million 
Monte Carlo moves with the accumulation of data at each temperature. Changes in the 
/^-parameter were made forty per cent of the time and the two-dimensional pressure was 
arbitrarily set to 1 atomic unit of pressure. Simulations were performed on 4x4, 8x8, 16x16 
and 32x32 lattices with periodic boundary conditions included. 

Shown in Figure is the constant pressure heat capacity of the compressible Ising model 
for 4x4, 8x8, 16x16 and 32x32 lattices from top to bottom in the figure. The peak at 240K 
in the heat capacity occurs at the same temperature as rapid changes in the magnetization 
of the model. As the number of spins in the simulation is increased the width of the peak 



narrows as expected. Given in Figure [14] is the isothermal compressibility as a function of 
temperature. Again from top to bottom is k for a 4x4, 8x8, 16x16 and 32x32 lattice. For 
the 4x4 lattice, no peak appears in the compressibility. As the sample size is increased, a 
peak develops in the compressibility at the transition temperature. Evidently, the existence 
of a peak in the compressibility is more sensitive to finite size effects than the heat capacity. 
Although not shown here, the isobaric coefficient of thermal expansion shows a peak at 
the transition temperature for all lattice sizes. These Ising results lend support to our 
assumption that the lack of peak in k in the isothermal-isobaric simulations of ammonium 
chloride is a result of finite size effects. 



19 



VI. SUMMARY AND DISCUSSION 



In this work we have applied Monte Carlo methods in the canonical and isothermal- 
isobaric ensembles to calculate the thermodynamic properties of NH 4 C1( S ) modeled by simple 
point charge pair interaction potentials. In both ensembles the model potential predicts an 
order to disorder transition associated with the rotational orientations of the ammonium ions 
in the lattice. In the isothermal-isobaric ensemble at 1 atmosphere pressure, the transition 
temperature is about 250K, as determined from the rapid growth of the susceptibility. Al- 
though we have made no effort to determine the transition temperature with any precision, 
the agreement between the approximately determined temperature and the experimental 
result (243K) is satisfying. 

The central sample cell used in the simulations contained a maximum of 27 ion pairs. 
Extensions to larger central cell sizes were inhibited by the large portion of the computer 
time needed to evaluate the Ewald corrections. Apparently this size limitation produced 
unreliable values for the isothermal compressibility, a quantity that compressible Ising cal- 
culations imply is sensitive to the size of the central cell. 

There are several outstanding questions about the behavior of ammonium chloride that 
require further attention. We found no evidence in our simulations for a rapid change in the 
molar volume at the transition temperature. Since the transition is known to be first order 
at atmospheric pressure jl9| , it is worth examining some of the physical effects not included 
in the simulations. A true discontinuity in the lattice constant would require an infinite 
central simulation cell. The unseen rapid change in volume may be another example of a 
finite size effect. However, in the simulations presented here we have included only those low 
frequency vibrational modes where the overall crystal symmetry is unaltered. Calculations 
that include more general variations in the locations of the ionic mass centers would clearly 
be of interest. Inclusion of the internal vibrations on the ammonium ions can be expected 
to be of less importance owing to the high frequencies of those motions. 

Another approximation in the results has been the application of classical mechanics. 
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Classical calculations can be interpreted as the infinite mass limit of a corresponding quan- 
tum calculation. If quantum effects were not important, our simulations would be equally 
applicable to the phase diagram of ND4C1( S ). Experimentally, a significant shift of the 
tricritical point to lower pressures is known to occur when the ammonium ions are deuter- 
ated |T%|Jim . It is possible that the classical system at atmospheric pressure has a continuous 
disordering transition, and the location (or existence) of a tricritical point is a consequence 
of quantum effects. A quantum path integral study of ammonium chloride within the same 
model potential would clearly be of interest. 
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TABLES 

TABLE I. The parameters used in the model potential 



Pair Aij C 6 D 12 

H-H a 1.0162 1.9950 2.9973 

N-N 6 104.74 1.5611 25.393 

Cl-Cl 6 125.55 1.7489 113.68 

H-N a 10.318 1.7780 8.7229 

H-C1 C 10.033 43884.0 

N-Cl a 114.22 1.6550 53.736 

a Combining rules 
b Reference [[32 1 



c Reference |}7 

d Units of energy in Hartree and units of length in Bohr 
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Figure drawn using XMol, version 1.3.1, Minnesota Supercomputer Center, Inc., Min- 
neapolis MN, 1993. 
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FIGURE CAPTIONS 



1. The two "spin" orientations of an ammonium ion in a cage of 8 chloride ions. Each 
orientation represents a potential minimum with the hydrogen atoms on the ammo- 
nium ion pointing to the shaded chloride ions. The two orientations are related to 
one another by rotating the ammonium ion by an angle of 7r/2 about the z-axis. The 
arrows between the hydrogen atoms and the chloride ions are included for clarity. |49 



2. The potential energy U in kJ mol _1 of an oriented ammonium chloride lattice with 27 
ion pairs in the central cell and Ewald corrections included. The energy is given as 
a function of the rotation angle <ft f° r a central ammonium ion as it is rotated out of 
orientation with the remaining ammonium ions in the lattice (see Figure 1). 

3. The difference in energy AU in kJ mol -1 between an ammonium chloride lattice with 
with all ammonium ions oriented in the same direction except the central ion and an 
ammonium chloride lattice with all ammonium ions oriented in the same direction 
except for the central ion and another ion a distance R in A from the central ion. The 
rapid approach to an asymptote is a consequence of the octapole-octapole interaction 
and lends support to the representation of the lattice by a nearest neighbor Ising 
model. 

4. The average potential energy in kJ mol _1 of NH4CI using 8 ion pairs in the central 
simulation cell as a function of temperature. Both curves were generated from an 
initial configuration with the ammonium ions in a random rotational configuration. 
In the upper curve ordinary Metropolis moves were used whereas in the lower curve a 
magnified step size was included ten percent of the time. The quasiergodicity problems 
evident in the upper curve at the lower temperatures are removed by the "magwalking" 
strategy. The error bars are smaller than the plotted points. 

5. An 8 ion pair representation of the ammonium chloride lattice. The representation 
is a particular configuration taken from a canonical simulation at 300K. The arrows 
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attached to each ammonium ion represent the algebraic sign of Mj. The up arrows 
are on ammonium ions with positive Mj and the down arrows are on ammonium ions 
with negative Mj. The arrows represent a mapping of the ammonium ion orientation 



onto a spin variable. E 



6. The constant volume heat capacity per ion pair in units of ks of ammonium chloride as 
a function of temperature. The data in the upper curve were obtained from canonical 
simulations with 8 ion pairs in the central cell, and the data in the lower curve were 
obtained from canonical simulations with 27 ion pairs in the central cell. By increasing 
the size of the central simulation cell the heat capacity narrows with a peak maximum 
at 210K. The error bars are at the double standard deviation level. 

7. The square of the magnetization as a function of temperature for a 27 ion pair rep- 
resentation of the ammonium chloride lattice in the canonical ensemble. The points 
represented by diamonds were computed directly using Eq.([L3|). The points repre- 
sented by triangles were obtained by mapping each ammonium ion orientation onto a 
spin variable as discussed in the text. The magnetization changes rapidly at the same 
temperature that the heat capacity shows a maximum in Figure |6|. The error bars in 
this figure are smaller than the representations of the plotted points. 

8. The susceptibility as a function of temperature for a 27 ion pair representation of the 
ammonium chloride lattice in the canonical ensemble. The error bars are at the double 
standard deviation level. The maximum in the peak of the susceptibility matches the 
maximum in the heat capacity and the region of rapid change in the magnetization. 

9. The square of the magnetization as a function of temperature for a 27 ion pair rep- 
resentation of the ammonium chloride lattice in the isothermal-isobaric ensemble at 
a pressure of 1 atmosphere. The points represented by diamonds were computed di- 
rectly using Eq.(|T^). The points represented by triangles were obtained by mapping 
each ammonium ion orientation onto a spin variable as discussed in the text. The 
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magnetization changes rapidly at the same temperature that the heat capacity and 
isobaric coefficient of thermal expansion have maxima. The error bars in this figure 
are smaller than the representations of the plotted points. 

10. The susceptibility as a function of temperature for a 27 ion pair representation of the 
ammonium chloride lattice in the isothermal-isobaric ensemble. The error bars are at 
the double standard deviation level. The maximum in the peak of the susceptibility 
identifies the onset of the order to disorder transition. 

11. The lattice parameter ao in A as a function of temperature for the 27 ion pair repre- 
sentation of the ammonium chloride crystal at 1 atmosphere pressure calculated in the 
isothermal-isobaric ensemble. Although no discontinuous change in the lattice param- 
eter is seen, a slope change is apparent at the transition temperature. The inflection 
matches the location of the peak in the isobaric coefficient of thermal expansion. The 
error bars are at the double standard deviation level. 

12. The constant pressure heat capacity C p per ion pair in units of ks, the isobaric 
coefficient of thermal expansion a in K _1 and the isothermal compressibility k in 
atmospheres" 1 as a function of temperature for the 27 ion pair representation of the 
ammonium chloride crystal at 1 atmosphere pressure in the isothermal-isobaric en- 
semble. The maxima in C p and a match the temperature at which there are rapid 
changes in the order parameter. The lack of a maximum in k is likely a result of finite 
size effects (see text). The error bars are at the double standard deviation level. 

13. The constant pressure heat capacity for the compressible Ising model. The four panels 
from top to bottom are results for 4x4, 8x8 16x16 and 32x32 lattices with periodic 
boundary conditions. The error bars are at the double standard deviation level. 

14. The isothermal compressibility in inverse atomic units of pressure for the compressible 
Ising model. The four panels from top to bottom are results for 4x4, 8x8 16x16 and 
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32x32 lattices with periodic boundary conditions. The error bars are at the double 
standard deviation level. 
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